Zenith Abacus Estimation
Template for RMD from mosaic; most code from MCF, some from JS
June 8, 2015
Load libraries
Read in the data first, then combine Woodcock and Wiat and z-score the Standardized Test composite
d <- read.csv("/Users/jsulliv1/Downloads/jesstimation-master/zenith all data.csv")
d$condition <- factor(d$condition, levels=c("control","abacus"))
d$standardized <- rowMeans(d[,c("wiat","woodcock")],na.rm=TRUE) Descriptives and preliminaries
Summaries.
d %>%
select(standardized, ans, deviance, linr2, ordinality, year) %>%
group_by(year) %>%
do(tidy(summary(.))) %>%
data.frame## year Var1 Var2 Freq
## 1 0 standardized Min. :0.04651
## 2 0 standardized 1st Qu.:0.15302
## 3 0 standardized Median :0.19302
## 4 0 standardized Mean :0.19596
## 5 0 standardized 3rd Qu.:0.23791
## 6 0 standardized Max. :0.34279
## 7 0 standardized <NA>
## 8 0 ans Min. :0.1000
## 9 0 ans 1st Qu.:0.2325
## 10 0 ans Median :0.3600
## 11 0 ans Mean :0.3657
## 12 0 ans 3rd Qu.:0.4600
## 13 0 ans Max. :0.7800
## 14 0 ans NA's :25
## 15 0 deviance Min. : NA
## 16 0 deviance 1st Qu.: NA
## 17 0 deviance Median : NA
## 18 0 deviance Mean :NaN
## 19 0 deviance 3rd Qu.: NA
## 20 0 deviance Max. : NA
## 21 0 deviance NA's :187
## 22 0 linr2 Min. : NA
## 23 0 linr2 1st Qu.: NA
## 24 0 linr2 Median : NA
## 25 0 linr2 Mean :NaN
## 26 0 linr2 3rd Qu.: NA
## 27 0 linr2 Max. : NA
## 28 0 linr2 NA's :187
## 29 0 ordinality Min. : NA
## 30 0 ordinality 1st Qu.: NA
## 31 0 ordinality Median : NA
## 32 0 ordinality Mean :NaN
## 33 0 ordinality 3rd Qu.: NA
## 34 0 ordinality Max. : NA
## 35 0 ordinality NA's :187
## 36 0 year Min. :0
## 37 0 year 1st Qu.:0
## 38 0 year Median :0
## 39 0 year Mean :0
## 40 0 year 3rd Qu.:0
## 41 0 year Max. :0
## 42 0 year <NA>
## 43 1 standardized Min. :0.07814
## 44 1 standardized 1st Qu.:0.26279
## 45 1 standardized Median :0.31116
## 46 1 standardized Mean :0.31115
## 47 1 standardized 3rd Qu.:0.36023
## 48 1 standardized Max. :0.54093
## 49 1 standardized <NA>
## 50 1 ans Min. :0.0200
## 51 1 ans 1st Qu.:0.1200
## 52 1 ans Median :0.1500
## 53 1 ans Mean :0.1773
## 54 1 ans 3rd Qu.:0.2000
## 55 1 ans Max. :0.7800
## 56 1 ans NA's :6
## 57 1 deviance Min. :0.230
## 58 1 deviance 1st Qu.:0.480
## 59 1 deviance Median :0.605
## 60 1 deviance Mean :0.708
## 61 1 deviance 3rd Qu.:0.880
## 62 1 deviance Max. :2.400
## 63 1 deviance NA's :3
## 64 1 linr2 Min. :0.0100
## 65 1 linr2 1st Qu.:0.2300
## 66 1 linr2 Median :0.3450
## 67 1 linr2 Mean :0.3663
## 68 1 linr2 3rd Qu.:0.5100
## 69 1 linr2 Max. :0.9100
## 70 1 linr2 NA's :3
## 71 1 ordinality Min. :0.5900
## 72 1 ordinality 1st Qu.:0.7575
## 73 1 ordinality Median :0.8000
## 74 1 ordinality Mean :0.7977
## 75 1 ordinality 3rd Qu.:0.8400
## 76 1 ordinality Max. :1.0000
## 77 1 ordinality NA's :3
## 78 1 year Min. :1
## 79 1 year 1st Qu.:1
## 80 1 year Median :1
## 81 1 year Mean :1
## 82 1 year 3rd Qu.:1
## 83 1 year Max. :1
## 84 1 year <NA>
## 85 2 standardized Min. :0.1298
## 86 2 standardized 1st Qu.:0.3744
## 87 2 standardized Median :0.4251
## 88 2 standardized Mean :0.4258
## 89 2 standardized 3rd Qu.:0.4801
## 90 2 standardized Max. :0.6358
## 91 2 standardized NA's :1
## 92 2 ans Min. :0.0100
## 93 2 ans 1st Qu.:0.1100
## 94 2 ans Median :0.1300
## 95 2 ans Mean :0.1462
## 96 2 ans 3rd Qu.:0.1700
## 97 2 ans Max. :0.5700
## 98 2 ans NA's :2
## 99 2 deviance Min. :0.2800
## 100 2 deviance 1st Qu.:0.4500
## 101 2 deviance Median :0.6100
## 102 2 deviance Mean :0.7127
## 103 2 deviance 3rd Qu.:0.8850
## 104 2 deviance Max. :3.0700
## 105 2 deviance <NA>
## 106 2 linr2 Min. :0.0600
## 107 2 linr2 1st Qu.:0.2050
## 108 2 linr2 Median :0.3200
## 109 2 linr2 Mean :0.3518
## 110 2 linr2 3rd Qu.:0.4700
## 111 2 linr2 Max. :0.8000
## 112 2 linr2 <NA>
## 113 2 ordinality Min. :0.6200
## 114 2 ordinality 1st Qu.:0.7500
## 115 2 ordinality Median :0.8000
## 116 2 ordinality Mean :0.7936
## 117 2 ordinality 3rd Qu.:0.8400
## 118 2 ordinality Max. :0.9700
## 119 2 ordinality <NA>
## 120 2 year Min. :2
## 121 2 year 1st Qu.:2
## 122 2 year Median :2
## 123 2 year Mean :2
## 124 2 year 3rd Qu.:2
## 125 2 year Max. :2
## 126 2 year <NA>
## 127 3 standardized Min. :0.1614
## 128 3 standardized 1st Qu.:0.4609
## 129 3 standardized Median :0.5442
## 130 3 standardized Mean :0.5381
## 131 3 standardized 3rd Qu.:0.6216
## 132 3 standardized Max. :0.7740
## 133 3 standardized <NA>
## 134 3 ans Min. :0.0600
## 135 3 ans 1st Qu.:0.1000
## 136 3 ans Median :0.1200
## 137 3 ans Mean :0.1384
## 138 3 ans 3rd Qu.:0.1600
## 139 3 ans Max. :0.5900
## 140 3 ans NA's :2
## 141 3 deviance Min. :0.2300
## 142 3 deviance 1st Qu.:0.4300
## 143 3 deviance Median :0.6250
## 144 3 deviance Mean :0.6984
## 145 3 deviance 3rd Qu.:0.9000
## 146 3 deviance Max. :2.0200
## 147 3 deviance NA's :1
## 148 3 linr2 Min. :0.0400
## 149 3 linr2 1st Qu.:0.1925
## 150 3 linr2 Median :0.3100
## 151 3 linr2 Mean :0.3560
## 152 3 linr2 3rd Qu.:0.4975
## 153 3 linr2 Max. :0.8200
## 154 3 linr2 NA's :1
## 155 3 ordinality Min. :0.5700
## 156 3 ordinality 1st Qu.:0.7525
## 157 3 ordinality Median :0.8000
## 158 3 ordinality Mean :0.7982
## 159 3 ordinality 3rd Qu.:0.8500
## 160 3 ordinality Max. :0.9500
## 161 3 ordinality NA's :1
## 162 3 year Min. :3
## 163 3 year 1st Qu.:3
## 164 3 year Median :3
## 165 3 year Mean :3
## 166 3 year 3rd Qu.:3
## 167 3 year Max. :3
## 168 3 year <NA>Reliability for Estimation tasks
Deviance.
data.frame(year=factor(c("1-2","2-3")),
corrs=c(cor.test(d$deviance[d$year==1],
d$deviance[d$year==2])$estimate,
cor.test(d$deviance[d$year==2],
d$deviance[d$year==3])$estimate))## year corrs
## 1 1-2 0.2775883
## 2 2-3 0.5415512data.frame(year=factor(c("1-2","2-3")),
pvals=c(cor.test(d$deviance[d$year==1],
d$deviance[d$year==2])$p.value,
cor.test(d$deviance[d$year==2],
d$deviance[d$year==3])$p.value))## year pvals
## 1 1-2 1.361945e-04
## 2 2-3 1.332268e-15Linear \(r^2\).
data.frame(year=factor(c("1-2","2-3")),
corrs=c(cor.test(d$linr2[d$year==1],d$linr2[d$year==2])$estimate,
cor.test(d$linr2[d$year==2],d$linr2[d$year==3])$estimate))## year corrs
## 1 1-2 0.1907387
## 2 2-3 0.3790020data.frame(year=factor(c("1-2","2-3")),
pvals=c(cor.test(d$linr2[d$year==1],d$linr2[d$year==2])$p.value,
cor.test(d$linr2[d$year==2],d$linr2[d$year==3])$p.value))## year pvals
## 1 1-2 9.499658e-03
## 2 2-3 9.585397e-08Ordinality.
data.frame(year=factor(c("1-2","2-3")),
corrs=c(cor.test(d$ordinality[d$year==1],
d$ordinality[d$year==2])$estimate,
cor.test(d$ordinality[d$year==2],
d$ordinality[d$year==3])$estimate))## year corrs
## 1 1-2 0.2130936
## 2 2-3 0.3334274data.frame(year=factor(c("1-2","2-3")),
pvals=c(cor.test(d$ordinality[d$year==1],
d$ordinality[d$year==2])$p.value,
cor.test(d$ordinality[d$year==2],
d$ordinality[d$year==3])$p.value))## year pvals
## 1 1-2 3.681035e-03
## 2 2-3 3.311927e-06ANS.
data.frame(year=factor(c("0-1","1-2","2-3")),
corrs=c(cor.test(d$ans[d$year==0],
d$ans[d$year==1])$estimate,
cor.test(d$ans[d$year==1],
d$ans[d$year==2])$estimate,
cor.test(d$ans[d$year==2],
d$ans[d$year==3])$estimate))## year corrs
## 1 0-1 0.2318435
## 2 1-2 0.1907536
## 3 2-3 0.4443862data.frame(year=factor(c("0-1","1-2","2-3")),
pvals=c(cor.test(d$ans[d$year==0],
d$ans[d$year==1])$p.value,
cor.test(d$ans[d$year==1],
d$ans[d$year==2])$p.value,
cor.test(d$ans[d$year==2],
d$ans[d$year==3])$p.value))## year pvals
## 1 0-1 3.377277e-03
## 2 1-2 1.010514e-02
## 3 2-3 2.630023e-10We see that the correlations are pretty decent year-to-year for our three DVs.
Training effects
OK, now we test for effect of training for ANS. We can use a model b/c we have all 4 years:
w.interaction <- lmer(ans ~ year * condition + (subnum|year), data=d)
summary(w.interaction)## Linear mixed model fit by REML ['lmerMod']
## Formula: ans ~ year * condition + (subnum | year)
## Data: d
##
## REML criterion at convergence: -1205.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.7297 -0.4898 -0.1617 0.3022 6.0360
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## year (Intercept) 4.513e-03 6.718e-02
## subnum 1.367e-12 1.169e-06 -1.00
## Residual 1.022e-02 1.011e-01
## Number of obs: 713, groups: year, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.324419 0.056844 5.707
## year -0.074453 0.030366 -2.452
## conditionabacus -0.021833 0.013056 -1.672
## year:conditionabacus 0.006231 0.006859 0.908
##
## Correlation of Fixed Effects:
## (Intr) year cndtnb
## year -0.802
## conditinbcs -0.113 0.091
## yr:cndtnbcs 0.092 -0.110 -0.814wo.interaction<-lmer(ans ~ year + condition + (subnum|year), data=d)
anova(w.interaction, wo.interaction)## refitting model(s) with ML (instead of REML)## Data: d
## Models:
## wo.interaction: ans ~ year + condition + (subnum | year)
## w.interaction: ans ~ year * condition + (subnum | year)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## wo.interaction 7 -1217.2 -1185.2 615.62 -1231.2
## w.interaction 8 -1216.0 -1179.5 616.03 -1232.0 0.8209 1 0.3649Now, we do t-tests by task to find differences between control and abacus groups. We will correct for multiple comparisons post-hoc, correcting for the # of comparisons per DV.
tasks <- c("deviance","linr2","ans", "ordinality")
d %>%
filter(year > 0) %>%
gather(measure, value, deviance, linr2, ans, ordinality) %>%
group_by(year, measure) %>%
do(tidy(t.test(.$value[.$condition == "abacus"],
.$value[.$condition == "control"])))## Source: local data frame [12 x 10]
## Groups: year, measure
##
## year measure estimate estimate1 estimate2 statistic p.value
## 1 1 deviance -0.080161367 0.6653488 0.7455102 -1.6349326 0.10385750
## 2 1 linr2 0.015889891 0.3747674 0.3588776 0.5947372 0.55279248
## 3 1 ans -0.002817625 0.1758140 0.1786316 -0.1979361 0.84333389
## 4 1 ordinality 0.021198386 0.8089535 0.7877551 2.1472626 0.03312995
## 5 2 deviance -0.117626263 0.6504545 0.7680808 -2.2387863 0.02637566
## 6 2 linr2 0.010303030 0.3572727 0.3469697 0.3780748 0.70582175
## 7 2 ans -0.002969775 0.1446591 0.1476289 -0.3397382 0.73446692
## 8 2 ordinality 0.014166667 0.8011364 0.7869697 1.3926937 0.16539755
## 9 3 deviance -0.072769070 0.6596552 0.7324242 -1.4796581 0.14067668
## 10 3 linr2 0.067293626 0.3918391 0.3245455 2.3208574 0.02146247
## 11 3 ans -0.011578700 0.1322989 0.1438776 -1.1261363 0.26159050
## 12 3 ordinality 0.011323581 0.8042529 0.7929293 1.1191883 0.26453958
## Variables not shown: parameter (dbl), conf.low (dbl), conf.high (dbl)PLOTS
ANS - Without intervention split
qplot(ans, standardized, facets=~year,
data=d) +
geom_smooth(method="lm") +
ylab("Standardized Test Composite") +
xlab("Weber Fraction") 
And with:
qplot(ans, standardized, facets=~year,
col=condition,
data=d) +
geom_smooth(method="lm") +
ylab("Standardized Test Composite") +
xlab("Weber Fraction") 
Now do the same for deviance: ```{/r, warnings=FALSE} qplot(deviance, standardized, facets=~year, data=d) + geom_smooth(method=“lm”) + ylab(“Standardized Test Composite”) + xlab(“PAE”)
qplot(deviance, standardized, facets=~year, col=condition, data=d) + geom_smooth(method=“lm”) + ylab(“Standardized Test Composite”) + xlab(“PAE”)
```
and for linear \(r^2\):
qplot(linr2, standardized, facets=~year,
data=d) +
geom_smooth(method="lm") +
ylab("Standardized Test Composite") +
xlab("Linear r^2") 
qplot(linr2, standardized, facets=~year,
col=condition,
data=d) +
geom_smooth(method="lm") +
ylab("Standardized Test Composite") +
xlab("Linear r^2") 
and for ordinality:
qplot(ordinality, standardized, facets=~year,
data=d) +
geom_smooth(method="lm") +
ylab("Standardized Test Composite") +
xlab("Ordinality") 
qplot(ordinality, standardized, facets=~year,
col=condition,
data=d) +
geom_smooth(method="lm") +
ylab("Standardized Test Composite") +
xlab("Ordinality") 
Models
Standardize predictors.
# this is insane, but PCA will choke if the scaled variables have included in them the scaling attributes (and this is what R does by default).
sscale <- function (x) {as.numeric(scale(x))}
md <- d %>%
gather(measure, value, deviance, linr2, ans, ordinality) %>%
group_by(year, measure) %>%
mutate(standardized.scale = sscale(standardized),
value.scale = sscale(value),
wiat.scale = sscale(wiat),
woodcock.scale = sscale(woodcock),
math.scale = sscale(math),
placeval.scale = sscale(placeval),
arith.scale = sscale(arith),
mental.rot.scale = sscale(mental.rot),
verbalwm.scale = sscale(verbalwm),
spatialwm.scale = sscale(spatialwm),
ravens.scale = sscale(ravens)) Now fit models. This uses the broom package - really very nifty.
NOTE, these are scaled values predicting scaled values.
std.baseline.models <- md %>%
group_by(year, measure) %>%
filter(!(year == 0 & measure != "ans")) %>%
do(tidy(lm(standardized.scale ~ value.scale, data = .))) %>%
filter(term != "(Intercept)") %>%
data.frame
std.baseline.models## year measure term estimate std.error statistic
## 1 0 ans value.scale -0.2861136935 0.07585437 -3.77188164
## 2 1 deviance value.scale -0.0736166477 0.07429885 -0.99081808
## 3 1 linr2 value.scale 0.1031406510 0.07410564 1.39180566
## 4 1 ans value.scale -0.1210668281 0.07315981 -1.65482699
## 5 1 ordinality value.scale 0.0008627001 0.07449894 0.01158003
## 6 2 deviance value.scale -0.1447405155 0.07301173 -1.98242838
## 7 2 linr2 value.scale 0.1910128279 0.07337888 2.60310359
## 8 2 ans value.scale -0.1049879875 0.07382586 -1.42210317
## 9 2 ordinality value.scale 0.2081562771 0.07257958 2.86797306
## 10 3 deviance value.scale -0.2283394329 0.07187531 -3.17688289
## 11 3 linr2 value.scale 0.2773777821 0.07093151 3.91050137
## 12 3 ans value.scale -0.2259644860 0.07165180 -3.15364723
## 13 3 ordinality value.scale 0.2275793284 0.07188841 3.16573055
## p.value
## 1 0.0002276970
## 2 0.3230898640
## 3 0.1656798389
## 4 0.0997119177
## 5 0.9907733616
## 6 0.0489197083
## 7 0.0099918231
## 8 0.1567077999
## 9 0.0046136812
## 10 0.0017460976
## 11 0.0001294263
## 12 0.0018849401
## 13 0.0018107944Check year 3 correlations
cor.test(d$standardized[d$year==3],d$linr2[d$year==3])##
## Pearson's product-moment correlation
##
## data: d$standardized[d$year == 3] and d$linr2[d$year == 3]
## t = 3.9105, df = 184, p-value = 0.0001294
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1386516 0.4047527
## sample estimates:
## cor
## 0.2770049cor.test(d$standardized[d$year==3],d$ans[d$year==3])##
## Pearson's product-moment correlation
##
## data: d$standardized[d$year == 3] and d$ans[d$year == 3]
## t = -3.1536, df = 183, p-value = 0.001885
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.35952906 -0.08557076
## sample estimates:
## cor
## -0.2270366cor.test(d$standardized[d$year==3],d$deviance[d$year==3])##
## Pearson's product-moment correlation
##
## data: d$standardized[d$year == 3] and d$deviance[d$year == 3]
## t = -3.1769, df = 184, p-value = 0.001746
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.36009736 -0.08700778
## sample estimates:
## cor
## -0.2280325cor.test(d$standardized[d$year==3],d$ordinality[d$year==3])##
## Pearson's product-moment correlation
##
## data: d$standardized[d$year == 3] and d$ordinality[d$year == 3]
## t = 3.1657, df = 184, p-value = 0.001811
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.08621322 0.35940039
## sample estimates:
## cor
## 0.2272734And here are the full models.
std.models <- md %>%
group_by(year, measure) %>%
filter(!(year == 0 & measure != "ans")) %>%
do(tidy(lm(standardized.scale ~ value.scale +
mental.rot.scale +
verbalwm.scale +
spatialwm.scale +
ravens.scale +
age + condition, data = .)))
options(dplyr.width = Inf)
print.data.frame(std.models)## year measure term estimate std.error statistic
## 1 0 ans (Intercept) -0.7856538941 0.97930903 -0.802253295
## 2 0 ans value.scale -0.2382604904 0.07869646 -3.027588537
## 3 0 ans mental.rot.scale -0.1506849858 0.08279761 -1.819919451
## 4 0 ans verbalwm.scale 0.1827323959 0.09210359 1.983987871
## 5 0 ans spatialwm.scale 0.2298199171 0.09478945 2.424530612
## 6 0 ans ravens.scale -0.0032556306 0.08233557 -0.039540997
## 7 0 ans age 0.1262516923 0.14629239 0.863009287
## 8 0 ans conditionabacus 0.0305031720 0.15180823 0.200932275
## 9 1 deviance (Intercept) 0.0704332170 0.87808443 0.080212352
## 10 1 deviance value.scale 0.0069503802 0.07071760 0.098283598
## 11 1 deviance mental.rot.scale 0.1684011132 0.08202233 2.053112944
## 12 1 deviance verbalwm.scale 0.1168479389 0.07678948 1.521665910
## 13 1 deviance spatialwm.scale 0.0875094913 0.07672455 1.140566990
## 14 1 deviance ravens.scale 0.2263509344 0.08458085 2.676148647
## 15 1 deviance age -0.0304643051 0.13145994 -0.231738320
## 16 1 deviance conditionabacus 0.2868962144 0.14741978 1.946117463
## 17 1 linr2 (Intercept) 0.0660351616 0.87843933 0.075173275
## 18 1 linr2 value.scale 0.0089746826 0.07048610 0.127325562
## 19 1 linr2 mental.rot.scale 0.1669835035 0.08228118 2.029425182
## 20 1 linr2 verbalwm.scale 0.1149181433 0.07730003 1.486650723
## 21 1 linr2 spatialwm.scale 0.0858767514 0.07573664 1.133886506
## 22 1 linr2 ravens.scale 0.2274777970 0.08449940 2.692064207
## 23 1 linr2 age -0.0297437056 0.13150432 -0.226180453
## 24 1 linr2 conditionabacus 0.2862859782 0.14722525 1.944544074
## 25 1 ans (Intercept) -0.1421182708 0.92169939 -0.154191564
## 26 1 ans value.scale -0.0848020275 0.06879832 -1.232617643
## 27 1 ans mental.rot.scale 0.1392803795 0.08123220 1.714595685
## 28 1 ans verbalwm.scale 0.1313498236 0.07610603 1.725879371
## 29 1 ans spatialwm.scale 0.0851861084 0.07508451 1.134536400
## 30 1 ans ravens.scale 0.2132227330 0.08329833 2.559748108
## 31 1 ans age 0.0095505003 0.13858947 0.068912163
## 32 1 ans conditionabacus 0.2263160447 0.14656650 1.544118489
## 33 1 ordinality (Intercept) 0.0876727683 0.87616117 0.100064659
## 34 1 ordinality value.scale -0.0646030826 0.07164768 -0.901677256
## 35 1 ordinality mental.rot.scale 0.1689033636 0.08172022 2.066849190
## 36 1 ordinality verbalwm.scale 0.1263959819 0.07720328 1.637184133
## 37 1 ordinality spatialwm.scale 0.0991810438 0.07684121 1.290727203
## 38 1 ordinality ravens.scale 0.2179444882 0.08476923 2.571033043
## 39 1 ordinality age -0.0335287511 0.13117888 -0.255595660
## 40 1 ordinality conditionabacus 0.2931081377 0.14707403 1.992929304
## 41 2 deviance (Intercept) -0.5018657569 0.88948097 -0.564223155
## 42 2 deviance value.scale -0.0511800660 0.07373884 -0.694071980
## 43 2 deviance mental.rot.scale 0.1230791522 0.08065492 1.525996802
## 44 2 deviance verbalwm.scale 0.1281293547 0.07541566 1.698975429
## 45 2 deviance spatialwm.scale -0.0340836956 0.07332083 -0.464856916
## 46 2 deviance ravens.scale 0.2693148811 0.08271429 3.255965458
## 47 2 deviance age 0.0580492047 0.13358344 0.434553904
## 48 2 deviance conditionabacus 0.2639209643 0.14321187 1.842870793
## 49 2 linr2 (Intercept) -0.3876979388 0.87850525 -0.441315448
## 50 2 linr2 value.scale 0.1380805185 0.07299928 1.891532724
## 51 2 linr2 mental.rot.scale 0.1124738978 0.08013608 1.403536268
## 52 2 linr2 verbalwm.scale 0.1099306123 0.07521069 1.461635399
## 53 2 linr2 spatialwm.scale -0.0389742974 0.07268825 -0.536184256
## 54 2 linr2 ravens.scale 0.2806527199 0.08098731 3.465391137
## 55 2 linr2 age 0.0396760723 0.13179446 0.301045078
## 56 2 linr2 conditionabacus 0.2764312569 0.14042647 1.968512480
## 57 2 ans (Intercept) -0.3072567066 0.94018041 -0.326806113
## 58 2 ans value.scale -0.0301034269 0.07696909 -0.391110614
## 59 2 ans mental.rot.scale 0.1146072990 0.08415605 1.361842738
## 60 2 ans verbalwm.scale 0.1365010632 0.07658785 1.782280867
## 61 2 ans spatialwm.scale -0.0410087042 0.07553182 -0.542932823
## 62 2 ans ravens.scale 0.2833881112 0.08283879 3.420958934
## 63 2 ans age 0.0270128435 0.14179116 0.190511472
## 64 2 ans conditionabacus 0.2836986921 0.14350583 1.976914100
## 65 2 ordinality (Intercept) -0.7205916314 0.88609373 -0.813222806
## 66 2 ordinality value.scale 0.1512482795 0.07301202 2.071553260
## 67 2 ordinality mental.rot.scale 0.1140304473 0.07989318 1.427286430
## 68 2 ordinality verbalwm.scale 0.1137262272 0.07457634 1.524963982
## 69 2 ordinality spatialwm.scale -0.0273514290 0.07247060 -0.377414128
## 70 2 ordinality ravens.scale 0.2668443782 0.08099826 3.294445631
## 71 2 ordinality age 0.0925382873 0.13313316 0.695080665
## 72 2 ordinality conditionabacus 0.2464247627 0.14098063 1.747933438
## 73 3 deviance (Intercept) -0.2560927282 0.80749465 -0.317144799
## 74 3 deviance value.scale -0.0357769329 0.06927171 -0.516472504
## 75 3 deviance mental.rot.scale 0.2709432040 0.07600035 3.565025514
## 76 3 deviance verbalwm.scale 0.1099748891 0.06626194 1.659699161
## 77 3 deviance spatialwm.scale 0.0994894786 0.07112983 1.398702634
## 78 3 deviance ravens.scale 0.2433194324 0.07590442 3.205603074
## 79 3 deviance age 0.0256576208 0.12083417 0.212337463
## 80 3 deviance conditionabacus 0.1907479142 0.13077011 1.458650670
## 81 3 linr2 (Intercept) -0.0655470230 0.80968285 -0.080953948
## 82 3 linr2 value.scale 0.1128325528 0.06946138 1.624392521
## 83 3 linr2 mental.rot.scale 0.2427626639 0.07753767 3.130899542
## 84 3 linr2 verbalwm.scale 0.0959844742 0.06643793 1.444724029
## 85 3 linr2 spatialwm.scale 0.1108330877 0.06911768 1.603541794
## 86 3 linr2 ravens.scale 0.2449441035 0.07497269 3.267110975
## 87 3 linr2 age -0.0008532564 0.12090824 -0.007057057
## 88 3 linr2 conditionabacus 0.1614555244 0.13088925 1.233527735
## 89 3 ans (Intercept) -0.3544991433 0.83241480 -0.425868383
## 90 3 ans value.scale -0.0629947326 0.06825250 -0.922965986
## 91 3 ans mental.rot.scale 0.2684758022 0.07572026 3.545627172
## 92 3 ans verbalwm.scale 0.1039107589 0.06691078 1.552974742
## 93 3 ans spatialwm.scale 0.0879826488 0.07247507 1.213971244
## 94 3 ans ravens.scale 0.2384733210 0.07604512 3.135944900
## 95 3 ans age 0.0412869395 0.12490245 0.330553476
## 96 3 ans conditionabacus 0.1844731761 0.13059066 1.412606203
## 97 3 ordinality (Intercept) -0.2044912280 0.80854880 -0.252911424
## 98 3 ordinality value.scale 0.0528422831 0.06870094 0.769163866
## 99 3 ordinality mental.rot.scale 0.2691781146 0.07574292 3.553838830
## 100 3 ordinality verbalwm.scale 0.1073684321 0.06635353 1.618126885
## 101 3 ordinality spatialwm.scale 0.0977280214 0.07056627 1.384911340
## 102 3 ordinality ravens.scale 0.2416496514 0.07578199 3.188747579
## 103 3 ordinality age 0.0180278107 0.12089009 0.149125624
## 104 3 ordinality conditionabacus 0.1896118295 0.13017170 1.456628716
## p.value
## 1 0.4237958964
## 2 0.0029461438
## 3 0.0709542685
## 4 0.0492556845
## 5 0.0166316159
## 6 0.9685166315
## 7 0.3896406260
## 8 0.8410495331
## 9 0.9361644093
## 10 0.9218249689
## 11 0.0416214241
## 12 0.1299840292
## 13 0.2556832131
## 14 0.0081894521
## 15 0.8170249786
## 16 0.0533190830
## 17 0.9401668218
## 18 0.8988359129
## 19 0.0440026059
## 20 0.1389926135
## 21 0.2584672279
## 22 0.0078244039
## 23 0.8213376636
## 24 0.0535099277
## 25 0.8776482769
## 26 0.2194825325
## 27 0.0883083896
## 28 0.0862525322
## 29 0.2582251573
## 30 0.0113779623
## 31 0.9451434874
## 32 0.1244870632
## 33 0.9204130166
## 34 0.3685269175
## 35 0.0402916937
## 36 0.1034749939
## 37 0.1985823287
## 38 0.0110117305
## 39 0.7985776068
## 40 0.0478980056
## 41 0.5733459095
## 42 0.4885844316
## 43 0.1288689680
## 44 0.0911530623
## 45 0.6426287760
## 46 0.0013640044
## 47 0.6644373417
## 48 0.0670893756
## 49 0.6595454225
## 50 0.0602538121
## 51 0.1622809782
## 52 0.1456871485
## 53 0.5925321054
## 54 0.0006702930
## 55 0.7637480085
## 56 0.0506350122
## 57 0.7442213401
## 58 0.6962110432
## 59 0.1750708691
## 60 0.0765095661
## 61 0.5878956429
## 62 0.0007835058
## 63 0.8491383711
## 64 0.0496877453
## 65 0.4172278081
## 66 0.0398175172
## 67 0.1553315807
## 68 0.1291262918
## 69 0.7063366207
## 70 0.0012000617
## 71 0.4879537205
## 72 0.0822809528
## 73 0.7515161450
## 74 0.6061839930
## 75 0.0004705116
## 76 0.0987866058
## 77 0.1636923296
## 78 0.0016053453
## 79 0.8320937454
## 80 0.1464744448
## 81 0.9355721535
## 82 0.1061126096
## 83 0.0020463495
## 84 0.1503441964
## 85 0.1106392139
## 86 0.0013103222
## 87 0.9943774607
## 88 0.2190522339
## 89 0.6707361801
## 90 0.3573178456
## 91 0.0005046884
## 92 0.1222672477
## 93 0.2264229451
## 94 0.0020151761
## 95 0.7413840986
## 96 0.1595784209
## 97 0.8006369287
## 98 0.4428448800
## 99 0.0004895490
## 100 0.1074570113
## 101 0.1678629307
## 102 0.0016963512
## 103 0.8816282024
## 104 0.1470314608gelb<-as.numeric(std.models$p.value)
small<-subset(std.models, gelb<.05)
print.data.frame(small)## year measure term estimate std.error statistic
## 1 0 ans value.scale -0.2382605 0.07869646 -3.027589
## 2 0 ans verbalwm.scale 0.1827324 0.09210359 1.983988
## 3 0 ans spatialwm.scale 0.2298199 0.09478945 2.424531
## 4 1 deviance mental.rot.scale 0.1684011 0.08202233 2.053113
## 5 1 deviance ravens.scale 0.2263509 0.08458085 2.676149
## 6 1 linr2 mental.rot.scale 0.1669835 0.08228118 2.029425
## 7 1 linr2 ravens.scale 0.2274778 0.08449940 2.692064
## 8 1 ans ravens.scale 0.2132227 0.08329833 2.559748
## 9 1 ordinality mental.rot.scale 0.1689034 0.08172022 2.066849
## 10 1 ordinality ravens.scale 0.2179445 0.08476923 2.571033
## 11 1 ordinality conditionabacus 0.2931081 0.14707403 1.992929
## 12 2 deviance ravens.scale 0.2693149 0.08271429 3.255965
## 13 2 linr2 ravens.scale 0.2806527 0.08098731 3.465391
## 14 2 ans ravens.scale 0.2833881 0.08283879 3.420959
## 15 2 ans conditionabacus 0.2836987 0.14350583 1.976914
## 16 2 ordinality value.scale 0.1512483 0.07301202 2.071553
## 17 2 ordinality ravens.scale 0.2668444 0.08099826 3.294446
## 18 3 deviance mental.rot.scale 0.2709432 0.07600035 3.565026
## 19 3 deviance ravens.scale 0.2433194 0.07590442 3.205603
## 20 3 linr2 mental.rot.scale 0.2427627 0.07753767 3.130900
## 21 3 linr2 ravens.scale 0.2449441 0.07497269 3.267111
## 22 3 ans mental.rot.scale 0.2684758 0.07572026 3.545627
## 23 3 ans ravens.scale 0.2384733 0.07604512 3.135945
## 24 3 ordinality mental.rot.scale 0.2691781 0.07574292 3.553839
## 25 3 ordinality ravens.scale 0.2416497 0.07578199 3.188748
## p.value
## 1 0.0029461438
## 2 0.0492556845
## 3 0.0166316159
## 4 0.0416214241
## 5 0.0081894521
## 6 0.0440026059
## 7 0.0078244039
## 8 0.0113779623
## 9 0.0402916937
## 10 0.0110117305
## 11 0.0478980056
## 12 0.0013640044
## 13 0.0006702930
## 14 0.0007835058
## 15 0.0496877453
## 16 0.0398175172
## 17 0.0012000617
## 18 0.0004705116
## 19 0.0016053453
## 20 0.0020463495
## 21 0.0013103222
## 22 0.0005046884
## 23 0.0020151761
## 24 0.0004895490
## 25 0.0016963512And make pretty plot.
std.models$term <- factor(std.models$term,
levels = c("(Intercept)",
"value.scale", "mental.rot.scale",
"spatialwm.scale", "verbalwm.scale",
"ravens.scale","age", "conditionabacus"),
labels = c("Intercept",
"Predictor", "Mental Rotation",
"Spatial WM", "Verbal WM",
"Raven's", "Age", "Intervention"))
std.models$measure <- factor(std.models$measure,
levels = c("ans","deviance","linr2","ordinality"),
labels = c("ANS","PAE","Linear r^2",
"Ordinality"))
qplot(term, estimate,
fill = term,
ymin = estimate - std.error, ymax = estimate + std.error,
geom = c("bar", "linerange"), stat = "identity",
facets = measure ~ year,
data=filter(std.models, term != "Intercept")) +
geom_hline(yintercept=0, lty=2) +
xlab("Predictor") +
ylab("Standardized Beta Weight") +
scale_fill_discrete(name="Predictor") +
theme(axis.text.x = element_text(angle = 90, vjust=.5, hjust = 1)) 
Add PCA to the dataset. Note that PC1 is flipped just to make life easier.
# need this to get around standard evals
pc1 <- function(x,y,z,m,n) {
sscale(as.numeric(prcomp(~x + y + z + m + n)$x[,1]))
}
# need to filter to complete cases
pmd <- md %>%
filter(complete.cases(wiat.scale, woodcock.scale, math.scale,
placeval.scale, arith.scale)) %>%
mutate(pc1 = -pc1(wiat.scale, woodcock.scale, math.scale,
placeval.scale, arith.scale))PCA Models and plot.
pca.baseline.models <- pmd %>%
group_by(year, measure) %>%
filter(!(year == 0 & measure != "ans")) %>%
do(tidy(lm(pc1 ~ value.scale, data = .))) %>%
filter(term != "(Intercept)") %>%
data.frame
pca.baseline.models## year measure term estimate std.error statistic
## 1 0 ans value.scale 0.2960511 0.07929057 3.733749
## 2 1 deviance value.scale -0.1676604 0.07346686 -2.282123
## 3 1 linr2 value.scale 0.1581167 0.07354619 2.149897
## 4 1 ans value.scale -0.1294271 0.07288281 -1.775825
## 5 1 ordinality value.scale 0.1054681 0.07430609 1.419374
## 6 2 deviance value.scale -0.1873279 0.07311760 -2.562009
## 7 2 linr2 value.scale 0.2497652 0.07302858 3.420102
## 8 2 ans value.scale -0.1474893 0.07411458 -1.990018
## 9 2 ordinality value.scale 0.1914047 0.07492045 2.554772
## 10 3 deviance value.scale -0.2476530 0.07146468 -3.465390
## 11 3 linr2 value.scale 0.3378703 0.06971814 4.846232
## 12 3 ans value.scale -0.2617409 0.07125388 -3.673356
## 13 3 ordinality value.scale 0.2735804 0.07165864 3.817828
## p.value
## 1 2.794896e-04
## 2 2.364572e-02
## 3 3.289062e-02
## 4 7.747032e-02
## 5 1.575100e-01
## 6 1.122584e-02
## 7 7.744574e-04
## 8 4.811134e-02
## 9 1.145338e-02
## 10 6.603674e-04
## 11 2.684857e-06
## 12 3.150773e-04
## 13 1.844388e-04pca.models <- pmd %>%
group_by(year, measure) %>%
filter(!(year == 0 & measure != "ans")) %>%
do(tidy(lm(pc1 ~ value.scale +
mental.rot.scale +
verbalwm.scale +
spatialwm.scale +
ravens.scale +
age + condition, data = .)))
options(dplyr.width = Inf)
print.data.frame(pca.models)## year measure term estimate std.error statistic
## 1 0 ans (Intercept) 1.497922211 0.97045769 1.54352140
## 2 0 ans value.scale 0.185976822 0.07924519 2.34685317
## 3 0 ans mental.rot.scale 0.072935906 0.08421307 0.86608773
## 4 0 ans verbalwm.scale -0.202937260 0.08948352 -2.26787308
## 5 0 ans spatialwm.scale -0.261960198 0.09449228 -2.77229205
## 6 0 ans ravens.scale -0.114758661 0.08801196 -1.30389846
## 7 0 ans age -0.229856115 0.14462246 -1.58935285
## 8 0 ans conditionabacus -0.093480549 0.15332473 -0.60968995
## 9 1 deviance (Intercept) -0.393256073 0.86642719 -0.45388243
## 10 1 deviance value.scale -0.076369306 0.06956242 -1.09785292
## 11 1 deviance mental.rot.scale 0.137902965 0.08062877 1.71034434
## 12 1 deviance verbalwm.scale 0.154817774 0.07547741 2.05118013
## 13 1 deviance spatialwm.scale 0.137621074 0.07537677 1.82577580
## 14 1 deviance ravens.scale 0.215687795 0.08310873 2.59524824
## 15 1 deviance age 0.041116779 0.12989342 0.31654243
## 16 1 deviance conditionabacus 0.260282871 0.14542286 1.78983458
## 17 1 linr2 (Intercept) -0.411056655 0.86683764 -0.47420259
## 18 1 linr2 value.scale 0.076027446 0.06928729 1.09727842
## 19 1 linr2 mental.rot.scale 0.134040458 0.08088218 1.65723111
## 20 1 linr2 verbalwm.scale 0.149613179 0.07597723 1.96918457
## 21 1 linr2 spatialwm.scale 0.147930968 0.07442020 1.98777981
## 22 1 linr2 ravens.scale 0.214289364 0.08303426 2.58073440
## 23 1 linr2 age 0.043047989 0.12994334 0.33128276
## 24 1 linr2 conditionabacus 0.269851367 0.14519564 1.85853635
## 25 1 ans (Intercept) -0.719905636 0.90360982 -0.79669966
## 26 1 ans value.scale -0.091798116 0.06736004 -1.36279788
## 27 1 ans mental.rot.scale 0.110370490 0.07934901 1.39094976
## 28 1 ans verbalwm.scale 0.179227511 0.07434468 2.41076429
## 29 1 ans spatialwm.scale 0.148767080 0.07333249 2.02866532
## 30 1 ans ravens.scale 0.193001655 0.08136940 2.37191943
## 31 1 ans age 0.099131013 0.13604215 0.72867867
## 32 1 ans conditionabacus 0.196803424 0.14371770 1.36937501
## 33 1 ordinality (Intercept) -0.387559399 0.86962288 -0.44566375
## 34 1 ordinality value.scale 0.023544324 0.07111527 0.33107270
## 35 1 ordinality mental.rot.scale 0.142111839 0.08079549 1.75890797
## 36 1 ordinality verbalwm.scale 0.157371241 0.07629260 2.06273283
## 37 1 ordinality spatialwm.scale 0.146513481 0.07591240 1.93003354
## 38 1 ordinality ravens.scale 0.212816809 0.08375460 2.54095678
## 39 1 ordinality age 0.039723398 0.13038597 0.30466006
## 40 1 ordinality conditionabacus 0.266280511 0.14594540 1.82452147
## 41 2 deviance (Intercept) -0.559636755 0.87668416 -0.63835618
## 42 2 deviance value.scale -0.072862999 0.07099878 -1.02625711
## 43 2 deviance mental.rot.scale 0.105192762 0.07799517 1.34870861
## 44 2 deviance verbalwm.scale 0.171976503 0.07220109 2.38191012
## 45 2 deviance spatialwm.scale -0.012544974 0.07086466 -0.17702723
## 46 2 deviance ravens.scale 0.301492650 0.07928570 3.80261067
## 47 2 deviance age 0.055304097 0.13162517 0.42016355
## 48 2 deviance conditionabacus 0.414240079 0.13779835 3.00613242
## 49 2 linr2 (Intercept) -0.440782477 0.85431608 -0.51594777
## 50 2 linr2 value.scale 0.191623271 0.06961279 2.75270222
## 51 2 linr2 mental.rot.scale 0.089100273 0.07673112 1.16120132
## 52 2 linr2 verbalwm.scale 0.145760406 0.07140294 2.04137808
## 53 2 linr2 spatialwm.scale -0.018563201 0.06949943 -0.26709862
## 54 2 linr2 ravens.scale 0.318710012 0.07676868 4.15156282
## 55 2 linr2 age 0.035902866 0.12811956 0.28022940
## 56 2 linr2 conditionabacus 0.433918051 0.13400545 3.23806259
## 57 2 ans (Intercept) -0.311710269 0.91191696 -0.34181870
## 58 2 ans value.scale -0.071477985 0.07407185 -0.96498183
## 59 2 ans mental.rot.scale 0.086794558 0.08049890 1.07820795
## 60 2 ans verbalwm.scale 0.174039099 0.07312008 2.38018203
## 61 2 ans spatialwm.scale -0.020299232 0.07245154 -0.28017670
## 62 2 ans ravens.scale 0.318265539 0.07901458 4.02793447
## 63 2 ans age 0.016343796 0.13727361 0.11906000
## 64 2 ans conditionabacus 0.434433459 0.13760369 3.15713528
## 65 2 ordinality (Intercept) -0.691272223 0.87443807 -0.79053308
## 66 2 ordinality value.scale 0.130566165 0.07138147 1.82913258
## 67 2 ordinality mental.rot.scale 0.101307039 0.07749056 1.30734685
## 68 2 ordinality verbalwm.scale 0.164190252 0.07171437 2.28950280
## 69 2 ordinality spatialwm.scale -0.003238722 0.07031302 -0.04606147
## 70 2 ordinality ravens.scale 0.305298433 0.07787643 3.92029326
## 71 2 ordinality age 0.075224262 0.13125771 0.57310357
## 72 2 ordinality conditionabacus 0.411803017 0.13615714 3.02446873
## 73 3 deviance (Intercept) 0.201030221 0.77841509 0.25825581
## 74 3 deviance value.scale -0.033882164 0.06661796 -0.50860401
## 75 3 deviance mental.rot.scale 0.265338664 0.07480796 3.54693091
## 76 3 deviance verbalwm.scale 0.093297204 0.06367089 1.46530389
## 77 3 deviance spatialwm.scale 0.133617372 0.06872634 1.94419458
## 78 3 deviance ravens.scale 0.288460659 0.07308947 3.94667858
## 79 3 deviance age -0.053056682 0.11630778 -0.45617481
## 80 3 deviance conditionabacus 0.328513763 0.12625532 2.60197950
## 81 3 linr2 (Intercept) 0.488534014 0.77079205 0.63380780
## 82 3 linr2 value.scale 0.174439128 0.06581358 2.65050345
## 83 3 linr2 mental.rot.scale 0.219014358 0.07509538 2.91648235
## 84 3 linr2 verbalwm.scale 0.070574026 0.06304830 1.11936438
## 85 3 linr2 spatialwm.scale 0.146887808 0.06588906 2.22932000
## 86 3 linr2 ravens.scale 0.288150617 0.07133822 4.03921794
## 87 3 linr2 age -0.092625073 0.11493704 -0.80587665
## 88 3 linr2 conditionabacus 0.278366611 0.12466335 2.23294668
## 89 3 ans (Intercept) 0.112100686 0.79960250 0.14019552
## 90 3 ans value.scale -0.091779712 0.06539806 -1.40340108
## 91 3 ans mental.rot.scale 0.260719124 0.07407022 3.51989114
## 92 3 ans verbalwm.scale 0.086171485 0.06410000 1.34432891
## 93 3 ans spatialwm.scale 0.118246779 0.06969060 1.69673923
## 94 3 ans ravens.scale 0.279641849 0.07308300 3.82635973
## 95 3 ans age -0.038693255 0.11982097 -0.32292558
## 96 3 ans conditionabacus 0.318335301 0.12554495 2.53562813
## 97 3 ordinality (Intercept) 0.283253764 0.77575484 0.36513309
## 98 3 ordinality value.scale 0.097249409 0.06583697 1.47712470
## 99 3 ordinality mental.rot.scale 0.258919954 0.07393976 3.50176886
## 100 3 ordinality verbalwm.scale 0.086623769 0.06349026 1.36436311
## 101 3 ordinality spatialwm.scale 0.124590058 0.06769418 1.84048404
## 102 3 ordinality ravens.scale 0.280743072 0.07272759 3.86020023
## 103 3 ordinality age -0.064727692 0.11582482 -0.55884127
## 104 3 ordinality conditionabacus 0.319673573 0.12493894 2.55863851
## p.value
## 1 1.254758e-01
## 2 2.065866e-02
## 3 3.882612e-01
## 4 2.522071e-02
## 5 6.503401e-03
## 6 1.948950e-01
## 7 1.147503e-01
## 8 5.432805e-01
## 9 6.505063e-01
## 10 2.738586e-01
## 11 8.907046e-02
## 12 4.182095e-02
## 13 6.968097e-02
## 14 1.029879e-02
## 15 7.519885e-01
## 16 7.530337e-02
## 17 6.359786e-01
## 18 2.741089e-01
## 19 9.936150e-02
## 20 5.059636e-02
## 21 4.848038e-02
## 22 1.072408e-02
## 23 7.408487e-01
## 24 6.486347e-02
## 25 4.267844e-01
## 26 1.748260e-01
## 27 1.661365e-01
## 28 1.703141e-02
## 29 4.411998e-02
## 30 1.886275e-02
## 31 4.672439e-01
## 32 1.727658e-01
## 33 6.564208e-01
## 34 7.410071e-01
## 35 8.043567e-02
## 36 4.069567e-02
## 37 5.530750e-02
## 38 1.197156e-02
## 39 7.610067e-01
## 40 6.987117e-02
## 41 5.241215e-01
## 42 3.062635e-01
## 43 1.792677e-01
## 44 1.835480e-02
## 45 8.597028e-01
## 46 2.008624e-04
## 47 6.749093e-01
## 48 3.057078e-03
## 49 6.065774e-01
## 50 6.568543e-03
## 51 2.472271e-01
## 52 4.279644e-02
## 53 7.897248e-01
## 54 5.267619e-05
## 55 7.796504e-01
## 56 1.453375e-03
## 57 7.329223e-01
## 58 3.359661e-01
## 59 2.825147e-01
## 60 1.844511e-02
## 61 7.796929e-01
## 62 8.567010e-05
## 63 9.053727e-01
## 64 1.894650e-03
## 65 4.303444e-01
## 66 6.917410e-02
## 67 1.929028e-01
## 68 2.330961e-02
## 69 9.633166e-01
## 70 1.291290e-04
## 71 5.673501e-01
## 72 2.886892e-03
## 73 7.965201e-01
## 74 6.116854e-01
## 75 5.030810e-04
## 76 1.446743e-01
## 77 5.351301e-02
## 78 1.155482e-04
## 79 6.488431e-01
## 80 1.008181e-02
## 81 5.270531e-01
## 82 8.791686e-03
## 83 4.015162e-03
## 84 2.645546e-01
## 85 2.709501e-02
## 86 8.086261e-05
## 87 4.214332e-01
## 88 2.685006e-02
## 89 8.886715e-01
## 90 1.623212e-01
## 91 5.541413e-04
## 92 1.806332e-01
## 93 9.157615e-02
## 94 1.824602e-04
## 95 7.471485e-01
## 96 1.212612e-02
## 97 7.154633e-01
## 98 1.414815e-01
## 99 5.896754e-04
## 100 1.742462e-01
## 101 6.743037e-02
## 102 1.604113e-04
## 103 5.770011e-01
## 104 1.137530e-02gelb2<-as.numeric(pca.models$p.value)
small2<-subset(pca.models, gelb2<.05)
print.data.frame(small2)## year measure term estimate std.error statistic
## 1 0 ans value.scale 0.1859768 0.07924519 2.346853
## 2 0 ans verbalwm.scale -0.2029373 0.08948352 -2.267873
## 3 0 ans spatialwm.scale -0.2619602 0.09449228 -2.772292
## 4 1 deviance verbalwm.scale 0.1548178 0.07547741 2.051180
## 5 1 deviance ravens.scale 0.2156878 0.08310873 2.595248
## 6 1 linr2 spatialwm.scale 0.1479310 0.07442020 1.987780
## 7 1 linr2 ravens.scale 0.2142894 0.08303426 2.580734
## 8 1 ans verbalwm.scale 0.1792275 0.07434468 2.410764
## 9 1 ans spatialwm.scale 0.1487671 0.07333249 2.028665
## 10 1 ans ravens.scale 0.1930017 0.08136940 2.371919
## 11 1 ordinality verbalwm.scale 0.1573712 0.07629260 2.062733
## 12 1 ordinality ravens.scale 0.2128168 0.08375460 2.540957
## 13 2 deviance verbalwm.scale 0.1719765 0.07220109 2.381910
## 14 2 deviance ravens.scale 0.3014926 0.07928570 3.802611
## 15 2 deviance conditionabacus 0.4142401 0.13779835 3.006132
## 16 2 linr2 value.scale 0.1916233 0.06961279 2.752702
## 17 2 linr2 verbalwm.scale 0.1457604 0.07140294 2.041378
## 18 2 linr2 ravens.scale 0.3187100 0.07676868 4.151563
## 19 2 linr2 conditionabacus 0.4339181 0.13400545 3.238063
## 20 2 ans verbalwm.scale 0.1740391 0.07312008 2.380182
## 21 2 ans ravens.scale 0.3182655 0.07901458 4.027934
## 22 2 ans conditionabacus 0.4344335 0.13760369 3.157135
## 23 2 ordinality verbalwm.scale 0.1641903 0.07171437 2.289503
## 24 2 ordinality ravens.scale 0.3052984 0.07787643 3.920293
## 25 2 ordinality conditionabacus 0.4118030 0.13615714 3.024469
## 26 3 deviance mental.rot.scale 0.2653387 0.07480796 3.546931
## 27 3 deviance ravens.scale 0.2884607 0.07308947 3.946679
## 28 3 deviance conditionabacus 0.3285138 0.12625532 2.601979
## 29 3 linr2 value.scale 0.1744391 0.06581358 2.650503
## 30 3 linr2 mental.rot.scale 0.2190144 0.07509538 2.916482
## 31 3 linr2 spatialwm.scale 0.1468878 0.06588906 2.229320
## 32 3 linr2 ravens.scale 0.2881506 0.07133822 4.039218
## 33 3 linr2 conditionabacus 0.2783666 0.12466335 2.232947
## 34 3 ans mental.rot.scale 0.2607191 0.07407022 3.519891
## 35 3 ans ravens.scale 0.2796418 0.07308300 3.826360
## 36 3 ans conditionabacus 0.3183353 0.12554495 2.535628
## 37 3 ordinality mental.rot.scale 0.2589200 0.07393976 3.501769
## 38 3 ordinality ravens.scale 0.2807431 0.07272759 3.860200
## 39 3 ordinality conditionabacus 0.3196736 0.12493894 2.558639
## p.value
## 1 2.065866e-02
## 2 2.522071e-02
## 3 6.503401e-03
## 4 4.182095e-02
## 5 1.029879e-02
## 6 4.848038e-02
## 7 1.072408e-02
## 8 1.703141e-02
## 9 4.411998e-02
## 10 1.886275e-02
## 11 4.069567e-02
## 12 1.197156e-02
## 13 1.835480e-02
## 14 2.008624e-04
## 15 3.057078e-03
## 16 6.568543e-03
## 17 4.279644e-02
## 18 5.267619e-05
## 19 1.453375e-03
## 20 1.844511e-02
## 21 8.567010e-05
## 22 1.894650e-03
## 23 2.330961e-02
## 24 1.291290e-04
## 25 2.886892e-03
## 26 5.030810e-04
## 27 1.155482e-04
## 28 1.008181e-02
## 29 8.791686e-03
## 30 4.015162e-03
## 31 2.709501e-02
## 32 8.086261e-05
## 33 2.685006e-02
## 34 5.541413e-04
## 35 1.824602e-04
## 36 1.212612e-02
## 37 5.896754e-04
## 38 1.604113e-04
## 39 1.137530e-02pca.models$term <- factor(pca.models$term,
levels = c("(Intercept)",
"value.scale", "mental.rot.scale",
"spatialwm.scale", "verbalwm.scale",
"ravens.scale","age", "conditionabacus"),
labels = c("Intercept",
"Predictor", "Mental Rotation",
"Spatial WM", "Verbal WM",
"Raven's", "Age", "Intervention"))
pca.models$measure <- factor(pca.models$measure,
levels = c("ans","deviance","linr2","ordinality"),
labels = c("ANS","PAE","Linear r^2",
"Ordinality"))
qplot(term, estimate,
fill = term,
ymin = estimate - std.error, ymax = estimate + std.error,
geom = c("bar", "linerange"), stat = "identity",
facets = measure ~ year,
data=filter(pca.models, term != "Intercept")) +
geom_hline(yintercept=0, lty=2) +
xlab("Predictor") +
ylab("Standardized Beta Weight") +
scale_fill_discrete(name="Predictor") +
theme(axis.text.x = element_text(angle = 90, vjust=.5, hjust = 1)) 
Additional code to double check for no glitches in above code
Now for the guts of one of the models, just an example from year 0, ANS. This uses ANOVA for model comparison.
## run models
# note, the filtering expression gets the model out of the same data frame that we used for the figure
# std.models <- md %>%
# group_by(year, measure) %>%
# filter(!(year == 0 & measure != "ans")) %>%
# do(tidy(lm(standardized.scale ~ value.scale +
# mental.rot.scale +
# verbalwm.scale +
# spatialwm.scale +
# ravens.scale +
# age + condition, data = .)))
model1 <- lm(standardized.scale ~ value.scale + mental.rot.scale + verbalwm.scale +
spatialwm.scale + ravens.scale + age + condition,
data = filter(md, year == 0, measure == "ans",
complete.cases(value.scale)))
model2 <- lm(standardized.scale ~ mental.rot.scale + verbalwm.scale +
spatialwm.scale + ravens.scale + age + condition,
data = filter(md, year == 0, measure == "ans",
complete.cases(value.scale)))
summary(model1)##
## Call:
## lm(formula = standardized.scale ~ value.scale + mental.rot.scale +
## verbalwm.scale + spatialwm.scale + ravens.scale + age + condition,
## data = filter(md, year == 0, measure == "ans", complete.cases(value.scale)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.13977 -0.63214 -0.01961 0.62713 2.29587
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.785654 0.979309 -0.802 0.42380
## value.scale -0.238260 0.078696 -3.028 0.00295 **
## mental.rot.scale -0.150685 0.082798 -1.820 0.07095 .
## verbalwm.scale 0.182732 0.092104 1.984 0.04926 *
## spatialwm.scale 0.229820 0.094789 2.425 0.01663 *
## ravens.scale -0.003256 0.082336 -0.040 0.96852
## age 0.126252 0.146292 0.863 0.38964
## conditionabacus 0.030503 0.151808 0.201 0.84105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9024 on 137 degrees of freedom
## (17 observations deleted due to missingness)
## Multiple R-squared: 0.2187, Adjusted R-squared: 0.1788
## F-statistic: 5.479 on 7 and 137 DF, p-value: 1.438e-05summary(model2)##
## Call:
## lm(formula = standardized.scale ~ mental.rot.scale + verbalwm.scale +
## spatialwm.scale + ravens.scale + age + condition, data = filter(md,
## year == 0, measure == "ans", complete.cases(value.scale)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.40108 -0.61155 -0.01836 0.60736 2.33124
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.178042 0.999003 -1.179 0.24034
## mental.rot.scale -0.140436 0.085141 -1.649 0.10133
## verbalwm.scale 0.156911 0.094382 1.663 0.09868 .
## spatialwm.scale 0.288327 0.095505 3.019 0.00302 **
## ravens.scale 0.009648 0.084623 0.114 0.90939
## age 0.182500 0.149340 1.222 0.22377
## conditionabacus 0.069035 0.155685 0.443 0.65815
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9287 on 138 degrees of freedom
## (17 observations deleted due to missingness)
## Multiple R-squared: 0.1665, Adjusted R-squared: 0.1302
## F-statistic: 4.593 on 6 and 138 DF, p-value: 0.0002778anova(model1, model2)## Analysis of Variance Table
##
## Model 1: standardized.scale ~ value.scale + mental.rot.scale + verbalwm.scale +
## spatialwm.scale + ravens.scale + age + condition
## Model 2: standardized.scale ~ mental.rot.scale + verbalwm.scale + spatialwm.scale +
## ravens.scale + age + condition
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 137 111.55
## 2 138 119.01 -1 -7.4635 9.1663 0.002946 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# also note that p value for the anova is the same as the coefficient p value in the models data frame
filter(std.models, year == 0, measure == "ANS")## Source: local data frame [8 x 7]
## Groups: year, measure
##
## year measure term estimate std.error statistic
## 1 0 ANS Intercept -0.785653894 0.97930903 -0.8022533
## 2 0 ANS Predictor -0.238260490 0.07869646 -3.0275885
## 3 0 ANS Mental Rotation -0.150684986 0.08279761 -1.8199195
## 4 0 ANS Verbal WM 0.182732396 0.09210359 1.9839879
## 5 0 ANS Spatial WM 0.229819917 0.09478945 2.4245306
## 6 0 ANS Raven's -0.003255631 0.08233557 -0.0395410
## 7 0 ANS Age 0.126251692 0.14629239 0.8630093
## 8 0 ANS Intervention 0.030503172 0.15180823 0.2009323
## p.value
## 1 0.423795896
## 2 0.002946144
## 3 0.070954268
## 4 0.049255685
## 5 0.016631616
## 6 0.968516632
## 7 0.389640626
## 8 0.841049533Year 0 Analyses by JS, to check betas
##
## Call:
## lm(formula = standardized.scale ~ ans.scale + mental.rot.scale +
## verbalwm.scale + spatialwm.scale + ravens.scale + age + condition,
## data = yc0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.15507 -0.63666 -0.01975 0.63162 2.31229
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.858627 0.987588 -0.869 0.38614
## ans.scale -0.239697 0.079171 -3.028 0.00295 **
## mental.rot.scale -0.146649 0.080580 -1.820 0.07095 .
## verbalwm.scale 0.178904 0.090174 1.984 0.04926 *
## spatialwm.scale 0.231096 0.095316 2.425 0.01663 *
## ravens.scale -0.003261 0.082482 -0.040 0.96852
## age 0.127154 0.147338 0.863 0.38964
## conditionabacus 0.030721 0.152894 0.201 0.84105
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9088 on 137 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.2187, Adjusted R-squared: 0.1788
## F-statistic: 5.479 on 7 and 137 DF, p-value: 1.438e-05##
## Call:
## lm(formula = standardized.scale ~ mental.rot.scale + verbalwm.scale +
## spatialwm.scale + ravens.scale + age + condition, data = yc0)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.4182 -0.6159 -0.0185 0.6117 2.3479
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.257841 1.007289 -1.249 0.21387
## mental.rot.scale -0.136675 0.082861 -1.649 0.10133
## verbalwm.scale 0.153623 0.092405 1.663 0.09868 .
## spatialwm.scale 0.289929 0.096035 3.019 0.00302 **
## ravens.scale 0.009665 0.084773 0.114 0.90939
## age 0.183804 0.150407 1.222 0.22377
## conditionabacus 0.069529 0.156799 0.443 0.65815
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9353 on 138 degrees of freedom
## (2 observations deleted due to missingness)
## Multiple R-squared: 0.1665, Adjusted R-squared: 0.1302
## F-statistic: 4.593 on 6 and 138 DF, p-value: 0.0002778## Analysis of Variance Table
##
## Model 1: standardized.scale ~ ans.scale + mental.rot.scale + verbalwm.scale +
## spatialwm.scale + ravens.scale + age + condition
## Model 2: standardized.scale ~ mental.rot.scale + verbalwm.scale + spatialwm.scale +
## ravens.scale + age + condition
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 137 113.15
## 2 138 120.72 -1 -7.5706 9.1663 0.002946 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Longitudial growth models
Let’s start by plotting some growth curves for standardized math measures.
qplot(year, standardized, facets=~subnum,
geom="line",
data=subset(d, subnum %in% unique(subnum)[1:30])) +
ylim(c(0,1)) +
ylab("Standardized Math Tests (WIAT & WJ)") +
xlab("Year")
Now add some models. Check whether linear or quadratic is more useful.
preds <- c("standardized","year","subnum")
d.cc <- d[complete.cases(d[,preds]),preds]
mod1 <- lmer(standardized ~ year + (year | subnum),
data=d.cc)
mod2 <- lmer(standardized ~ year + I(year^2) + (year + I(year^2) | subnum),
data=d.cc)
anova(mod1,mod2)## refitting model(s) with ML (instead of REML)## Data: d.cc
## Models:
## mod1: standardized ~ year + (year | subnum)
## mod2: standardized ~ year + I(year^2) + (year + I(year^2) | subnum)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## mod1 6 -1931.1 -1903.4 971.54 -1943.1
## mod2 10 -1928.2 -1882.0 974.10 -1948.2 5.1145 4 0.2758No gain for this DV (standardized tests) in adding a quadratic term. Let’s drop it then.
Here’s a basic model showing individaul variability plotted by the same model for all participants. It’s a basic linear model and it doesn’t do that badly.
The approach we’re going to follow is to try and see gains by adding predictors to this model to better capture individual variation. (Note that if we use individualized growth terms in a mixed model we already fit pretty much all the variance so there’s not much else we can do, so we are using regular linear models here).
preds <- c("standardized","year","subnum")
d.cc <- d[complete.cases(d[,preds]),preds]
mod <- lm(standardized ~ year,
data=d.cc)
d.cc$mod <- predict(mod)
qplot(year, standardized, facets=~subnum,
geom="point",
data=subset(d.cc, subnum %in% unique(subnum)[1:30])) +
geom_line(aes(x=year, y=mod),lty=2) +
ylim(c(0,1)) +
ylab("Standardized Math Tests (WIAT & WJ)") +
xlab("Year")
Let’s try adding concurrent ANS as a predictor. As you can see, it does very little here if we have a single coefficient.
(Note that there are some missing values for ANS from year 0, so that’s why curves are truncated).
preds <- c("standardized","year","subnum","ans")
d.cc <- d[complete.cases(d[,preds]),preds]
mod <- lm(standardized ~ year,
data=d.cc)
mod2 <- lm(standardized ~ year + ans,
data=d.cc)
d.cc$mod <- predict(mod)
d.cc$mod2 <- predict(mod2)
anova(mod,mod2)## Analysis of Variance Table
##
## Model 1: standardized ~ year
## Model 2: standardized ~ year + ans
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 710 5.2014
## 2 709 5.0825 1 0.11887 16.582 5.184e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1qplot(year, standardized, facets=~subnum,
geom="point",
data=subset(d.cc, subnum %in% unique(subnum)[1:30])) +
geom_line(aes(x=year, y=mod),lty=2) +
geom_line(aes(x=year, y=mod2),lty=1,col="red") +
ylim(c(0,1)) +
ylab("Standardized Math Tests (WIAT & WJ)") +
xlab("Year")
Multiple coefficients, one for each year, and you start to see some teensy-tiny shape differences between curves. (I added abacus condition here as well).
The dashed line is without ANS, the solid is with.
preds <- c("standardized","year","subnum","condition","ans")
d.cc <- d[complete.cases(d[,preds]),preds]
mod <- lm(standardized ~ factor(year) * condition,
data=d.cc)
mod2 <- lm(standardized ~ factor(year) * condition +
factor(year) * ans - ans,
data=d.cc)
d.cc$mod <- predict(mod)
d.cc$mod2 <- predict(mod2)
anova(mod,mod2)## Analysis of Variance Table
##
## Model 1: standardized ~ factor(year) * condition
## Model 2: standardized ~ factor(year) * condition + factor(year) * ans -
## ans
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 704 5.1389
## 2 700 4.9528 4 0.18608 6.5748 3.386e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1qplot(year, standardized, facets=~subnum,
geom="point",col=condition,
data=subset(d.cc, subnum %in% unique(subnum)[1:30])) +
geom_line(aes(x=year, y=mod),lty=2) +
geom_line(aes(x=year, y=mod2),lty=1) +
ylim(c(0,1)) +
ylab("Standardized Math Tests (WIAT & WJ)") +
xlab("Year")
Just for kicks, let’s add everything into the mix and see how we do.
We are doing a bunch better in making customized predictions. In terms of significant effects, we get:
- arithmetic (big predictor of - you guessed it - standardized arithmetic)
- abacus condition (whew)
- mental rotation
- raven’s
- verbal wm (but not spatial)
Nothing for ANS. Note that these predictions get noticeably worse when you leave out arithmetic scores. Domain knowledge is a big predictor.
preds <- c("standardized","year","subnum","condition","ans","mental.rot","spatialwm","verbalwm","ravens","age")
d.cc <- d[complete.cases(d[,preds]),preds]
mod <- lm(standardized ~ factor(year) * condition,
data=d.cc)
mod2 <- lm(standardized ~ factor(year) + condition +
age + mental.rot +
spatialwm + ravens + verbalwm + ans,
data=d.cc)
d.cc$mod <- predict(mod)
d.cc$mod2 <- predict(mod2)
summary(mod2)##
## Call:
## lm(formula = standardized ~ factor(year) + condition + age +
## mental.rot + spatialwm + ravens + verbalwm + ans, data = d.cc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.32374 -0.05260 -0.00326 0.05339 0.20364
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.076577 0.045052 1.700 0.089649 .
## factor(year)1 0.078369 0.010729 7.305 8.01e-13 ***
## factor(year)2 0.177043 0.011678 15.161 < 2e-16 ***
## factor(year)3 0.280678 0.012651 22.185 < 2e-16 ***
## conditionabacus 0.013745 0.006179 2.224 0.026461 *
## age 0.004044 0.005990 0.675 0.499858
## mental.rot 0.069274 0.019269 3.595 0.000348 ***
## spatialwm 0.005252 0.003459 1.518 0.129385
## ravens 0.084440 0.016781 5.032 6.27e-07 ***
## verbalwm 0.013164 0.003741 3.519 0.000463 ***
## ans -0.078625 0.031119 -2.527 0.011748 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07883 on 662 degrees of freedom
## Multiple R-squared: 0.7316, Adjusted R-squared: 0.7276
## F-statistic: 180.5 on 10 and 662 DF, p-value: < 2.2e-16qplot(year, standardized, facets=~subnum,
geom="point",col=condition,
data=subset(d.cc, subnum %in% unique(subnum)[1:30])) +
geom_line(aes(x=year, y=mod),lty=2) +
geom_line(aes(x=year, y=mod2),lty=1) +
ylim(c(0,1)) +
ylab("Standardized Math Tests (WIAT & WJ)") +
xlab("Year")
A final analysis: Let’s try this by lagging the predictors. We use the year before to predict the current year - so this is trying to use e.g. year 0’s ANS to see if it predicts growth over standard predictions in year 1.
Note that I added gender to the baseline predictor set so we will see only lag differences.
Define a function to do the lagging.
lag <- function(x, p) {
for (y in x$year) {
if (y == 0) {
eval(parse(text=paste("x$",p,".lag[x$year==0] <- 0",sep="")))
}
if (y != 0 & ((y-1) %in% x$year)) {
eval(parse(text=paste("x$",p,".lag[x$year==y] <- x$",
p,"[x$year==y-1]",sep="")))
} else {
eval(parse(text=paste("x$",p,".lag[x$year==y] <- NA",sep="")))
}
}
return(x)
}Now do the actual predictions. We are actually missing too many observations to really do this right without some serious imputation. We lose around a third of our data, as you can see from the plot.
Because of the missing data problem, do this with just ANS, to see if anything else is different. Results:
- ANS is a predictor
- lagged ANS isn’t
- None of this is visible in the plot because the differences are too small, so I left out the plot.